///<author>Arjen van Zanten</author>
///<datecreated>2012-03-27</datecreated>
///<summary>A static class that holds a number of Otsu functions.</summary>
///<review>Theodoor de Graaff
    ///<date>28-3-2012</date>
///</review>

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Drawing;
using System.Drawing.Imaging;
using KentekenHerkenning;

namespace KentekenHerkenning
{
    public class Otsu
    {
        /// <summary>
        /// Gets a threshold using a single pass of the standard Otsu's method.
        /// </summary>
        /// <param name="histArray">The input histogram array</param>
        /// <returns>The determined threshold</returns>
        public static int SinglePass(int[] histArray)
        {
            return GetThreshold(histArray, 0);
        }

        /// <summary>
        /// Gets a threshold using a double pass of the standard Otsu's method, where the second pass is applied to the part of the histogram below the threshold found in the first pass.
        /// </summary>
        /// <param name="histArray">The input histogram array</param>
        /// <returns>The determined threshold</returns>
        public static int DoublePassLow(int[] histArray)
        {
            int threshold = GetThreshold(histArray,0);

            int[] subHistArray = new int[threshold];
            Array.Copy(histArray, subHistArray, threshold);

            return GetThreshold(subHistArray, 0);
        }

        /// <summary>
        /// Gets a threshold using a double pass of the standard Otsu's method, of which the second pass is applied below the threshold found in the firrst pass.
        /// </summary>
        /// <param name="histArray">The input histogram array</param>
        /// <returns>The determined threshold</returns>
        public static int DoublePassHigh(int[] histArray)
        {
            int threshold = GetThreshold(histArray, 0);

            int[] subHistArray = new int[histArray.Length - threshold];
            Array.Copy(histArray, threshold, subHistArray, 0, histArray.Length - threshold);

            return GetThreshold(subHistArray, threshold);
        }

        /// <summary>
        /// Gets a threshold using a triple pass of the standard Otsu's method, of which the second pass is applied below the threshold found in the first pass, and the third pass applied between the first and second threshold.
        /// </summary>
        /// <param name="histArray"></param>
        /// <returns></returns>
        public static int TriplePassLowHigh(int[] histArray)
        {
            int threshold1 = GetThreshold(histArray, 0);
            int[] histArray1 = new int[threshold1];
            Array.Copy(histArray, histArray1, threshold1);
            int threshold2 = GetThreshold(histArray1, 0);

            int[] histArray2 = new int[threshold1 - threshold2];
            Array.Copy(histArray1, threshold2, histArray2, 0, threshold1 - threshold2);
            return GetThreshold(histArray2, threshold2);
        }

        /// <summary>
        /// Gets a threshold using a single pass of the standard Otsu's method, with an offset on the bin numbers.
        /// </summary>
        /// <param name="histArray">The input histogram array</param>
        /// <param name="offset">The offset</param>
        /// <returns>The determined threshold</returns>
        public static int GetThreshold(int[] histArray, int offset)
        {
            int[] cluster;
            double[] classVarianceArray = new double[histArray.Length]; // array that will hold all the within class viarances 
            int finalThreshold = 0;

            for (int threshold = 0; threshold < histArray.Length; threshold++)
            {
                cluster = new int[threshold];
                Array.Copy(histArray, 0, cluster, 0, threshold); // copy everything below the threshold to the cluster array
                classVarianceArray[threshold] += getWeightedVariance(cluster, (double)histArray.Sum(), offset); // add the weighted variance of the first cluster to the Within Class Variance total for this threshold.

                cluster = new int[histArray.Length - threshold];
                Array.Copy(histArray, threshold, cluster, 0, histArray.Length - threshold); // copy everything above the threshold to the cluster array
                classVarianceArray[threshold] += getWeightedVariance(cluster, (double)histArray.Sum(), offset + threshold); // add the weighted variance of the second cluster to the Within Class Variance total for this threshold.
            }

            // find the threshold with the smalles Within Class Variance, ignoring 0 or below.
            double min = -1;
            for (int i = 0; i < classVarianceArray.Length; i++)
            {
                if (classVarianceArray[i] > 0 && (min == -1 || classVarianceArray[i] < min))
                {
                    min = classVarianceArray[i];
                    finalThreshold = i;
                }
            }

            return finalThreshold + offset;
        }


        /// <summary>
        /// Gets the weighted variance for a given cluster, with a given histogram sum and bin offset.
        /// </summary>
        /// <param name="cluster">Input cluster array</param>
        /// <param name="histSum">Sum of the histogram</param>
        /// <param name="offset">Bin offset</param>
        /// <returns>Weighted variance</returns>
        private static double getWeightedVariance(int[] cluster, double histSum, int offset)
        {
            // if the cluster has 0 bins or all the bins are empty, return -1, so it can be ignored. 
            if (cluster.Length == 0 || cluster.Sum() == 0)
            {
                return -1;
            }

            // Calculate the mean
            double mean = 0; 
            for (int bin = 0; bin < cluster.Length; bin++)
            {
                mean += (bin + offset) * cluster[bin];
            }
            mean /= (double)cluster.Sum();

            // Calculate the variance
            double variance = 0;
            for (int bin = 0; bin < cluster.Length; bin++)
            {
                variance += Math.Pow((bin + offset - mean), 2) * cluster[bin];
            }
            variance /= histSum; // SHORTCUT: weight is actually not needed (surprisingly). 

            return variance;
        }
    }
}